#ifndef BIO_PH_INITIALIZER
#define BIO_PH_INITIALIZER

#include "kernel/chemsys.hpp"

//! Affine constraint on some chemical concentrations
class ChemIni
{
public:
	const Real value;
	
	//! a weight and a name for a species
	class Term
	{
	public:
		inline Term( const Real w, const string &n ) : weight( w ), name( n ) {}
		inline ~Term() throw() {}
		Term( const Term &t ) : weight( t.weight ), name( t.name ) {}
		
		const Real   weight;
		const string name;
		
	private:
		YOCTO_DISABLE_ASSIGN(Term);
	};
	
	ChemIni( const Library &lib, const Real V );
	~ChemIni() throw();
	ChemIni( const ChemIni &ci );
	
	void append( const Term &t );

	bool is_same_as( const ChemIni &other ) const throw();
	
private:
	vector<Term> terms_;
	YOCTO_DISABLE_ASSIGN(ChemIni);
	friend class Initializer;
public:
	const Library &library;
};


//! Supplemental conditions to solve the chemical system
class Initializer
{
public:
	explicit Initializer( const Library &library );
	virtual ~Initializer() throw();
	
	const size_t Nc; //!< #initializer
	const size_t M;  //!< #species
	const size_t N;  //!< #would be reactions = M-Nc
	
	void append( const ChemIni &ci );
	void load( lua_State *L, const string &ini_name); //!< load from a Lua VM
	void compile();
	void setup( ChemSys &chemsys, Real t );

private:
	vector<ChemIni>  I;     //!< all the constraints
	matrix<Real>     C;     //!< dI/dX, NcxM, Nc can be 0, precomputed at compile()
	//matrix<Real>     P;     //!< MxNc, PseudoInverse of C = C'*inv(C*C')
	matrix<Real>     Q;     //!< MxN,  in orthogonal of C
	//vector<Real>     V0;    //!< 1..Nc = I.value
	vector<Real>     X0;    //!< P * V0, default X value, which happens to have the minimum module
	linsys<Real>     lss;   //!< linear solver, at most MxM
	YOCTO_DISABLE_COPY_AND_ASSIGN(Initializer);
	
public:
	const Library &library;
	Real           C0; //!< default concentration for zero values of X0 (Y0=Q'*[conc_ini..conc_ini]') (1e-5)
	
};

#endif
